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ABSTRACT 

We calculate the energy spectra of cosmic ray (CR) protons and electrons at a plane shock with 
quasi-parallel magnetic fields, using time-dependent, diffusive shock acceleration (DSA) simulations, 
including energy losses via synchrotron emission and Inverse Compton (IC) scattering. A thermal 
leakage injection model and a Bohm type diffusion coefficient are adopted. The electron spectrum at 
the shock becomes steady after the DSA energy gains balance the synchrotron/IC losses, and it cuts 
off at the equilibrium momentum p^q- In the postshock region the cutoff momentum of the electron 
spectrum decreases with the distance from the shock due to the energy losses and the thickness of the 
spatial distribution of electrons scales as p~^ ■ Thus the slope of the downstream integrated spectrum 
steepens by one power of p for pbr < p < Pcq, where the break momentum decrease with the shock 
age as pbr oc In a CR modified shock, both the proton and electron spectrum exhibit a concave 
curvature and deviate from the canonical test-particle power-law, and the upstream integrated electron 
spectrum could dominate over the downstream integrated spectrum near the cutoff momentum. Thus 
the spectral shape near the cutoff of X-ray synchrotron emission could reveal a signature of nonlinear 
DSA. 
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1. INTRODUCTION 

Diffusive shock acceleration (DSA) is widely ac- 
cepted as the primary mechanism through which cos- 
mic rays (CRs) are produced in a variety of astrophys- 
ical environments (Blandford & Eichler 1987; Drury 
1983). Detailed nonlinear treatments of DSA predict 
that a small fraction of incoming thermal particles can 
be injected into the CR population, and accelerated to 
very high energies through their interactions with reso- 
nantly scattering Alfven waves in the converging flows 
across collisionless shocks (e.g. Berezhko, Ksenofontov 
& Volk 2002, 2009; Kang et al. 2002). 

Since DSA operates on relativistic particles with 
the same rigidity {R = pc/Ze) in exactly the same 
manner, both electrons and protons are expected to 
be accelerated at shocks. Indeed, multi-band obser- 
vations of nonthermal radio to 7-ray emissions from 
several supernova remnants have confirmed the accel- 
eration of CR electrons and protons up to 100 TeV (e.g. 
Volk, Berezhko, & Ksenofontov 2005; Parizot et al. 
2006; Abdo et al. 2010). On the other hand, post- 
shock thermal electrons need to be pre-accelerated be- 
fore they can be injected into Fermi process, since elec- 
trons have smaller gyroradius than protons by a fac- 
tor of (me/mp)^/^ (Reynods 2008). This results in 



a much smaller injection rate for electrons, i.e., the 
CR electron-to-proton ratio is estimated to be i^e/p ^ 
10~* — 10~^, which is not yet constrained accurately 
by plasma physics. Moreover, electrons lose energy pri- 
marily by synchrotron emission and Inverse Compton 
(IC) scattering, and so their energy spectrum has a 
cutoff much lower than that of the proton spectrum. 

In a strong, non-relativistic shock, DSA theory pre- 
dicts that the power-law distribution function for elec- 
trons, fe{p) oc with g « 4, which translates into the 
synchrotron flux, S„ oc with a = (g — 3)/2 « 0.5. 
For a unresolved, optically-thin source, we should ob- 
serve the volume integrated emission spectrum that 



steepens to v' 



-(g-2)/2 



OC V 



because the volume inte- 
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grated electron spectrum steepens to Fe{p) (x 
above a break momentum, pbr, due to synchrotron/IC 
losses. Such characteristics of nonthermal electrons and 
their synchrotron emission have been explored exten- 
sively in several previous studies (e.g. Webb et al. 1984; 
Heavens & Meisenheimer 1987; Berezhko, Ksenofontov 
& Volk 2002; Zirakashvili & Aharonian 2007). 

Webb et al. (1984) and Heavens & Meisenheimer 
(1987) considered steady-state solutions for the electron 
spectrum accelerated at a plane shock in the case of 
constant, momentum-independent diffusion coefHcient, 
while Zirakashvili & Aharonian (2007) generalize it to 
the case of momentum-dependent diffusion coefficient. 
In the present paper, we performed time-dependent nu- 
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merical simulations in which DSA of protons and elec- 
trons are followed along with synchrotron/IC losses for 
the shock parameters relevant for typical young super- 
nova remnants (SNRs). A thermal leakage injection 
model and a Bohm-type diffusion coefficient cx p) 

are adopted. Even though the CR proton spcctnmi ex- 
tends to ever higher momentum with time, the electron 
spectrum at the shock approaches to time-asymptotic 
states and can be compared with the steady-state so- 
lution of Zirakashvili & Aharonian (2007). We also 
demonstrate how the nonlinear feedback of the CR 
pressure affects the electron spectrum in a CR mod- 
ified shock. 

In the next section we describe the numerical method 
and model parameters. Some analytic estimates for 
the various features in the electron spectrum will be 
presented in Section 3. The simulation results will be 
discussed in Section 4, followed by a brief summary in 
Section 5. 

2. NUMERICAL METHOD 



2.1 CRASH code for DSA 

Here we consider the CR acceleration at quasi- 
parallel shocks where the magnetic field lines are par- 
allel to the shock normal. So we solve the stan- 
dard gasdynamic equations with CR proton pressure 
terms added in the conservative, Eulcrian formulation 
for one dimensional plane-parallel geometry. The ba- 
sic gasdynamic equations and details of the CRASH 
(Cosmic-Ray Amr SHock) code for one-dimensional 
plane-parallel geometry can be found in (Kang et al. 
2002). 

We solve the following diffusion-convection equa- 
tions for the pitch-angle-averaged phase space distri- 
bution function for CR protons, fp{x,p,t), and for CR 
electron, fe{x,p,t) (SkiUing 1975): 
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where ft = p'*/p, ft = p'^fe, y = ln(p). Here the 
particle momentum is expressed in units of nipC and 
so the spatial diffusion coefficient, k{x,p), has the 
same form for both protons and electrons. The cool- 
ing term b{p) = —dp/dt takes account for electron syn- 
chrotron/IC losses. 



The velocity represents the effective relative mo- 
tion of scattering centers with respect to the bulk fiow 
velocity, u. The mean wave speed is set to the Alfven 
speed, i.e., = va = Bj ^Amp in the upstream re- 
gion. This term reflects the fact that the scattering 
by Alfven waves tends to isotropizc the CR distribu- 
tion in the wave frame rather than the gas frame (Bell 
1978). In the postshock region, = is assumed, 
since the Alfvenic turbulence in that region is proba- 
bly relatively balanced (Jones 1993). This reduces the 
velocity difference between upstream and downstream 
scattering centers compared to the bulk flow, leading 
to less efficient DSA. Hereafter we use the subscripts 
'0', '1', and '2' to denote conditions far upstream of the 
shock, immediately upstream of the gas subshock and 
immediately downstream of the subshock, respectively. 

The dynamical effects of the CR proton pressure are 
included in the DSA simulations, while the CR elec- 
trons are treated as test-particles. In equations (1) 
and (2) the spatial advection terms arc solved by the 
wave-propagation algorithm, while the diffusion term 
is solved by the Crank-Nicholson scheme (Kang et al. 
2002). 

The radiative cooling term in equation (2) is treated 
by the operator splitting method, in which the follow- 
ing advection equation in the momentum y space is 
solved: 

where the advection speed is V = —[b{p)/p^] ■ p = 
—DB^p. Here the synchrotron cooling constant is de- 
fined as D = 4e'*/(9mgC^) (in cgs imits), and B^ = 
{B + B'^y/'^ is the "effective" magnetic field strength 
and B'^/S-K is the energy density of ambient radiation 
field (Longair 1994). The downward advection in mo- 
mentum introduces an additional Courant condition for 
time step, M < ^y/\V\. 

2.2 Model Parameters 

Here we consider a plane-parallel shock with Us = 
3000 km s~^, propagating into the upstream medium 
with the following properties: gas temperature, Tq = 
10^ K, sound speed, Cg = 150 km s^^ (sonic Mach 
number, M = Ug/cs = 20), and hydrogen number 
density, no = lcm~^. Because high- resolution X-ray 
observations of several young SNRs indicate the pres- 
ence of magnetic fields as strong as a few 100/zG down- 
stream of the shock [e.g., Parizot et aZ.2006), we adopt 
Bq = 30 for the upstream field strength, which is 
stronger than the mean ISM field of 5-8 jiG. The up- 
stream Alfven speed is va = Bo/y/Airpo = 65.5 km s~^ 
(Alfvenic Mach number, Ma = Us/va = 45.8) and 
so va/cs = 0.44. This set of parameters may repre- 
sent young supernova remnants in the hot interstel- 
lar medium with self-amplified magnetic fields (Kang 
2010). 

The physical quantities are normalized, both in the 
numerical code and in the plots below, by the following 
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constants: 



3.1 Basics of Diffusive Shock Acceleration 



Uo ^ 3000 km s-\ 
Ko = 1-04 X lO^^cm^s"^ 
To = Ko/uo = 3.47 X lO^^cm, 
to = I'^o/ul = 366years. 

Note that the calculations depend on the density po 

only through Alfvcn speed. So for the cases where 
Alfvenic drift is ignored, the density parameter uq is 
arbitrary. 

Since we consider only quasi-parallel shocks, the 
magnetic field is passive in the simulations. However, 

theoretical studies have shown that efficient magnetic 
field amplification via resonant and non-resonant wave- 
particle interactions is an integral part of DSA (Lucek 
& Bell 2000, Bell 2004), and the evolution of magnetic 
field strength is important in modeling the diffusion 
coefficient and synchrotron cooling. For simplicity, we 
assume the local magnetic field strength scales with 
gas density as B{x,t) = Bo[p{x,t)/po\, which repre- 
sents enhancement of resonant Alfven wave amplitudes 
through compression. 

We adopt a Bohm-type diffusion coefficicint that 
includes a weaker non-relativistic momentum depen- 
dence, 



P 



p{x) 



Po 



(4) 



where the coefficient k* 
lO^icm^s-i (So/30 AtG)-^ 



TOpc3/(3eBo) = 1-04 x 
The density dependence 
approximately accounts for the compressive amplifica- 
tion of Alfven waves. 

Also we adopt a phenomenological injection model, 
in which particles above a certain injection momentum 
Pinj cross the shock and get injected to the CR pop- 
ulation (Kang et al. 2002). The "effective" injection 
momentum can be approximated by 



Pinj w 1.17mpU2 (l + 



(5) 



where pih = yj2mpkBT2 is the thermal peak momen- 
tum of the immediate postshock gas with tempera- 
ture T2 and ks is the Boltzmann constant (Kang & 
Ryu 2010). One free parameter controls this function; 
es = Bo/B±, the ratio of the general magnetic field 
along the shock normal, Bq, to the amplitude of the 
postshock MHD wave turbulence, B±. Since we are 
not interested in the absolute amplitude of the elec- 
tron spectrum in this study, protons and electrons are 
injected in the same manner in the simulations, ignor- 
ing the fact that postshock thermal electrons have a 
much smaller gyroradius. In other words, we take ef- 
fectively the electron-to-proton ratio, i^e/p = 1- 



In the test-particle regime of DSA theory, the post- 
shock CR spectrum takes the power-law form of f{p) oc 
where the spectral index, q, depends on the veloc- 
ity jump across the shock: 



q = 



3(ui - va) 



VA-U2 



(6) 



where Ui = Us and U2 are the upstream and down- 
stream speed, respectively. Alfvenic drift of self- 
generated waves in the upstream region reduces the 
velocity jump that the particles experience across the 
shock, which in turn softens the CR spectrum beyond 
the canonical test-particle slope. 

The mean acceleration time for a particle to reach 
p from an injection momentum pinj in the test-particle 
limit of DSA theory is given by (Drury 1983) 



tacc(p) 



Ul — VA — U2 



Pinj 
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(7) 

Then the mean acceleration time scale is estimated by 



rririC 



(8) 



where F 3ct(2 - M^i)/(1 - MX^){a - 1 - aM^^), 
a = Ui/u2 = P2/P1 is the shock compression ra- 
tio, and Ma = ui/va is Alfven Mach number (Kang 
& Ryu 2010). In the limit of large Ma, the factor 
F « 6a/{(7 - 1) = 2(7, and F « 8 for cr = 4. So the 
maximum momentum accelerated by the shock age of t 
can be estimated aspmax/wpC « (uj./ FK*)t, which rep- 
resents the cutoff momentum in the proton spectrum. 
For typical supernova remnant shocks, it becomes 



Pn 



3.3x10^ ( 



t 



Bp 

3000 km s- 1 y V yrs y V 30 /uG 

(9) 



3.2 Radiative Cooling 

Relativistic electrons lose energy primarily through 
combined synchrotron and IC radiation with the cool- 
ing rate b{p). So the cooling time scale is 



P 



t.M = ^^={DBlp)-' 



1.33 X 10^ yrs 



D \ -2 / \ -1 

Be \ IP 



100 jiG J \mpC 



(10) 



This can be rewritten in term of the electron Lorentz 
factor, 7e = p/rUeC as 



ANALYTIC ESTIMATES 



irad(7e) = 2.45x10" yrs 
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The electron cutoff energy can be estimated from 
the equiUbrium condition that the momentum gain per 
cycle by DSA is equal to the synchrotron/IC loss per 
cycle, i.e., (Ap)DSA + (Ap)iad=0. Following the deriva- 
tion given by Webb et al. (1984), 



(Ap)dsa 

(Ap)rad = - 



Ap {ui - 1*2 - va) 
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where v is the particle velocity, and 



(12) 
(13) 

c,l - ^1 ^- Br 

and i?c,2 = Bi + B"^. The upstream and downstream 
residence time per cycle are A^i = 4ki/(ui — va)v 
and At2 = 4:K2/u2V, respectively. For ki — aK2 = 
(c^ /ieBi) ■ p, the equilibrium momentum for electrons 
is given in cgs units by 
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If we assume B2 — 
(7 = 4, this becomes 
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(15) 

which does not depend on the shock age t unlike the 
proton cutoff momentum Pmax- 

For p <^ pcq the DSA momentum gain is much larger 
than the radiative losses, so electrons and protons are 
accelerated in almost the same manner until the radia- 
tive loss balances the DSA gain. This occurs at the 
equilibrium time, teq = iacc(Poq), after which the elec- 
tron spectrum at the shock becomes steady and cuts 
off at Pcq. Using equations (8) and (15) with F = 8, it 
gives 



icq ~ (71.7 yrs) 
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Figure 1 shows the evolution of the proton and 
electron spectrum at the shock position from the test- 
particle DSA simulation (see next section for details). 
The electron cutoff momentum asymptotes to Peq after 
t ^ 0.5to ~ 183 yrs, a few times longer than teq, while 
the proton maximum momentum continues to increase 
with time as Pmax c>c t. For p < peq, the distribution 
function, gp{xs,p) and Qeixs^p) evolve in exactly the 
same way, as expected. 

Using equation (15), we can estimate the character- 
istic synchrotron frequency of electrons with p^q as 
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Fig. 1. — Time Evolution of the proton and electron 
distribution functions, gp{xs,p) (black dotted lines) and 
ge{xs,p) (red solid lines), at the shock position at t/to = 
0.1,0.2,0.3,0.5, 1., 3. from the test-particle DSA simula- 
tion for Us = 3000 km s"\ Bo = 30 //G. The Alfenic drift 
is not included {va = 0). 



which depends only on the shock speed. So the turn- 
over frequency of X-ray synchrotron emission from 
SNRs provides a measure of the shock speed, but not 
the magnetic field strength directly. 

Zirakashvili & Aharonian (2007) suggested that the 
electron spectrum at the shock position for a steady- 
state plane shock can be approximated by 

fe{xs,p) cc p'^[l + 0.523(p/po)'/']' exp( V/Po). 

(18) 

in the case of Bohm diffusion with B2 = VTIBq and 
(7 = 4. Their cutoff momentum po is given by 



/ B2 



1.56 X 10^ / Us 
g(l -hCT-i/2) V3OOO km s-iy VlOO fiG 



-1/2 



Po 

which is similar to, but not the same as, Pcq in equation 
(15), 

At the shock the electron spectrum cuts off at Peq, 
but in the postshock region electrons lose energy while 
being advected downstream, and so the electron cutoff 
momentum moves downward in the momentum space. 
The cutoff momentum, Pcut, at the distance d down- 
stream from the shock is given by 



Pcut(rf) 



Pcq 



1 + d/l 



(20) 



cq 



where l^q = U2 • irad(Poq) = U2[DBl.2Pcq) '■ In 
the limit of d 3> Icq, the cutoff momentum of the 
electron spectrum decreases away from the shock as 
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Fig. 2. — Particle spectra from the test-particle DSA simulation for Us = 3000 km s~^, Bo = 30 fiG. Top: Proton 
spectrum gp{x,p) (dotted lines) and electron spectrum ge{x,p) (solid lines) at 5 different locations downstream of the shock 
at t/to = 5 are shown. The curves for the electron spectrum are labeled with the distance downstream from the shock in 
the number of grid zones (-100, -50, -10, -5, 0). Note that the cooling length for electrons with Pcq, ^eq, corresponds to 6.15 
grid zones in the simulation. The dashed line shows the analytic form for the electron spectrum at the shock (Eq. [18]) 
suggested by Zirakashvili & Aharonian (2007). Bottom: Same as above except the spectra at two locations (-f 10 and +50 
grid zones) upstream of the shock are shown along with the spectra at the shock position. 



Pcut — Pcq • (Icq/d) oc d^^ . Thls Implles that the thick- 
ness of the spatial distribution of electrons with mo- 
mentum p decrease as Ax{p) — U2trad(p) (x p~^. 

At the shock age t ^ t^q, the electron spectrum at 
the furthest downstream point (dadv — U2 ■ t) has the 

cutoff at Pcut = Phl{t) « Pcq{lcq/U2t) W {DBl ^t)''^ . 

This gives 

^ « 1.34 X 10^ ( ( -^"l ' , (21) 
rupc \ 103 yrs ) \im ^iG ) ' ^ ' 

which depends only on the postshock magnetic field 
strength and the shock age, but not on the shock speed. 
If the electron distribution function at the shock has the 
form, fe(xs,p) oc for p < pcq, then the downstream 
integrated electron spectrum is given by 

/■° 

F,a{t,p)= /e(x,p,i)dxcxp-(«+i), (22) 



for pbr < p < Pcq, which was also shown in previ- 
ous studies (e.g. Heavens & Meisenheimer 1987; Zi- 
rakashvili & Aharonian 2007). 

4. DSA SIMULATION RESULTS 
4.1 Test-Particle Case 

In order to make a direct comparison between the 
analytic estimates and the DSA simulation results, we 
first consider a test-particle case, in which the dynam- 
ical feedback of the CR pressure is ignored. So the 
shock structure is determined by the canonical shock 
jump conditions. The injection parameter, eg = 0.2 is 
adopted and the Alfvenic drift is ignored {i.e., va — 0) 
for this case. 

Figure 2 shows the proton and electron spectra at 
different locations at t = 5to = 1.83 x 10^ yrs. It shows 
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Fig. 3. — Proton and electron distribution function, gp{x,p) and ge{x,p), in tlie pliase-space at t/to = 0.4, 1, 3 for tlie 
test-particle case siiown in Fig. 2. At eacii contour level, the value of g increases by a factor of 10. Note that the spatial 
span of the proton distribution shown here is [-4,-1-4], while that of the electron distribution is [-2,-f2]. 



that the cutoff momentum for the proton spectrum is 
Pmax/nipC sa 6 X 10^ (Eq. [9]), while that for the elec- 
tron spectrum at the shock is Peq/mpC « 2.4 x 10''^ (Eq. 
[15]). For comparison, we also show the analytic form 
in equation (18) suggested by Zirakashvili & Aharonian 
(2007), which agrees well with the simulated electron 
spectrum at the shock, especially near the exponential 
cutoff. 

The DSA operates efficiently only within a diffu- 
sion length, ^ k(p)/us, from the shock. Far behind 
the shock radiative cooling of electrons dominates over 
DSA, so the function ge{x,p) shifts downward in Inp 
space in a manner approximately described by the ad- 
vection equation in (3). Because the advection speed 
depends on electron momentum as V (x p, the shape 
of the cutoff steepens as we move downstream from 
the shock. So the electron spectrum well behind shock 
cannot be approximated by a simple analytic form 
such as a power-law with a Gaussian exponential cut- 
off. Heavens & Meisenheimer (1987) also observed the 
same trend in the electron spectrum behind a steady- 
state plane shock in the case of constant magnetic field 
strength and momentum-independent diffusion coeffi- 



cient {i.e., Bi — i?2 and ki — K2). According to 
equation (20) the cutoff momentum is estimated to be 
Pcut/mpC « 1.3 X 10^, 2.5 x 10^, 8.8 x 10^, 1.3 x 10^ at 
100, 50, 10, 5 grid zones, respectively, downstream from 
the shock. These estimates match well the simulated 
electron spectra, as can be seen in Figure 2. 

In upstream region, only the highest energy particles 
diffuse away from the shock, while low energy particles 
are swept downstream. So there the proton spectrum 
gp{p) has a peak at Pmax, while the electron spectrum 
ge{p) has a peak at pcq, as can be seen in the bottom 
panel of Figure 2. The upstream electron spectrum has 
the same exponential cutoff at p^q as the spectrum at 
the shock. Again the electron spectrum has the same 
shape as the proton spectrum for p < p^q. 

Figure 3 shows the evolution of the phase-space dis- 
tribution of gp{x,p) and ge{x,p) for the test-particle 
case. At low energy, p/nipC < 10"^, electron cooling is 
not important and so protons and electrons have the 
same distribution. At higher energy, electrons have 
much narrower distribution both in x and Inp space 
due to radiative cooling. It also demonstrates that the 
spatial thickness of the electron distribution narrows 
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Fig. 4. — Top: DSA simulation results from the test-particle case: Us = 3000 km s~^, Bo ~ 30 fiG, va = 0, es = 
0.2. Bottom: DSA simulation results from the CR modified case: Us ~ 3000 km s~^, Bo = 30 /iG, va = 65.5 km s~^, 
es = 0.25. Proton and electron distribution functions, gp{x,p) (black solid lines) and ge{x,p) (red lines), at the shock 
location (left panels) and at a upstream location (middle panels) at t/to = 1, 3, 5. The right panels show the downstream 

integrated electron spectrum, Ge,2 = fe{p)p'^dx (red dashed lines), the upstream integrated electron spectrum, Ge,i = 
/o^°° /e(p)p*da; (blue dotted lines) and the sum of the two curves (green solid lines) at t/to ~ 1, 3, 5. 



down at higher momentum as Ax (x p^^ in the post- 
shock region, and that the cutoff momentum Pcut{x) 
decreases as the downstream distance increases away 
from the shock. Although the postshock electron dis- 
tribution extends downstream as the flow advects, the 
upstream electron distribution does not change after 
t/to ^ 1 due to the balance between DSA and cooling. 
On the contrary, high energy protons diffuse much fur- 
ther both downstream and upstream as Pmax oc t in- 
creases with time. 

The top panels of Figure 4 show the evolution of 
gp{p) and geip) at the shock (right panel) and at a up- 
stream position (middle panel) atf = 3.7xl0^, 1.1 x 
10^, 1.8 X 10^ yrs. Because t t^^ ^ 71.7 yrs, so 
the electron spectrum has already reached the steady- 
state, as can be seen in the figure. In the right panel we 
also plot the downstream integrated electron spectrum, 

Ge,2 = /"^oo 5e(p)da;, thc upstream integrated electron 
spectrum, Ge^i = ge{p)dx. As discussed right af- 



ter equation (16), Ge,2 (dashed lines) becomes a broken 
power-law which steepens from p^'^ to ^"(9+1) above 
P > Phi (t) with an exponential cutoff at the higher mo- 
mentum peq ~ 2.4 X 10*. We also see that the brake mo- 
mentum decrease with time as pbr oc t^^. Note that the 
upstream integrated spectrum Ge,i (dotted lines) has 
reached the steady-state, while, for p < pbr, the am- 
plitude of the downstream spectrum Ge.2 increase lin- 
early with time as the flow advects downstream. Such 
a trend was seen in Figure 3 as well. Total volume 
integrated spectrum (solid lines) shows a small bump 
near the cutoff momentum due to the upstream contri- 
bution, which in turn will determine the exact shape of 
the cutoff of X-ray synchrotron emission. 

4.2 CR Modified Case 

With the injection parameter es = 0.25, the CR 
injection and acceleration is efficient enough to mod- 
ify significantly the shock structure. The CR injection 
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U3=3000km/s, B(, = 30/zG, modified 




Fig. 5. — Proton and electron distribution functions, gp(x,p) and ge{x,p), in the pliase-space at t/to = 0.4, 1, 3 for the 
CR modified case shown in the bottom panels of Fig. 4. At each contour level, the value of g increases by a factor of 10. 
Note that the spatial span of the proton distribution shown here is [-4, +4], while that of the electron distribution is [-2, +2]. 



fraction becomes ^ « 5 x 10^'*, and the postshock CR 
pressure is Pcr.2/ Poul ~ 0.29. So the postshock gas 
density, P2/P0 ~ 5.6, is larger than cr = 4 for the gasdy- 
namic shock. In such CR modified shocks, the pressure 
from CRs diffusing upstream compresses and deceler- 
ates the gas smoothly before it enters the subshock, 
creating a shock precursor (Kang & Jones 2007). 

With the assumed momentum-dependent diffusion, 
the particles of different momenta, p, experi- 
ence different compressions, depending on their diffu- 
sion length, ld{p) = k{p)/us. The particles just above 
Pinj sample mostly the compression across the sub- 
shock (cTg = pi/pq), while those near Pmax experience 
the total compression across the entire shock struc- 
ture {at = P2/ Po)- This leads to the particle distribu- 
tion function that behaves as f{p) oc p~3crs/(crs-i) 
p ^ pinj, but flattens gradually to f{p) oc p^^'^t/i<rt-i) 
toward p ~ Pmax (Kang et al. 2009). 

The bottom panels of Figure 4 show the evolution 
of gp{p), ge{p) at the shock position (left panel) and at 
a upstream location (middle panel) and the volume in- 
tegrated Ge{p) (right panel) as discussed before. The 
Alfvenic drift with va — 0.44cs {Ma = 45) are con- 



sidered. Because of the development of a smooth pre- 
cursor and the weaker subshock, both gp and g^ at 
the shock are softer than the test-particle power-law 
spectra at lower momenta p/{mpc) < 10^, while they 
are harder at higher momenta with a cutoff at Pmax 
(proton spectrum) or p^q (electron spectrum). Thus 
the CR spectrum exhibits the well-known concave cur- 
vature between the lowest and the highest momenta. 
Such concavity is reflected in the volume integrated 
spectrum as well, so Ge is no longer a simple broken 
power-law as in the test-particle case. 

With the greater velocity jump {at — 5.6), the ac- 
celeration is more efficient and so there are more high- 
est energy particles in the upstream region, compared 
to the test-particle case. As a result, the upstream 
integrated spectrum, Ge,i (dotted lines) has a more 
pronounced peak at peq, compared to the test-particle 
case. This introduces an additional curvature in the to- 
tal Ge spectrum. In fact, Ge,i dominates over Ge,2 near 
Pcq, so the upstream contribution should determine the 
spectral shape of X-ray synchrotron emission. Thus the 
spectral slope in radio and the detail shape in X-ray of 
the observed synchrotron flux can provide a measure of 
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nonlinear DSA feedback. 

Finally, Figure 5 shows the phase-space distribution 
of gp{x,p) and ge{x,p) for the CR modified model with 
high injection rate. Because the postshock magnetic 
field strength, B2 = (JtBo = 5.6i3o, is stronger, elec- 
trons cool down to lower energies, compared to the 
test-particle case shown in Figure 3. Enhanced cool- 
ing also reduces the thickness of the electron spatial 
distribution downstream of the shock. Again, we can 
see that at the highest energies of p/nipC > lO**, the 
upstream electron components is more important than 
the downstream component. 

5. SUMMARY 

Using the kinetic simulations of difFiisivc shock ac- 
celeration at a plane shock, we calculate the time- 
dependent evolution of the CR proton and electron 
spectra, including electronic synchrotron/IC energy 
losses. Both protons and electrons are injected at the 
shock via thermal leakage injection and accelerated by 
DSA, while electrons are treated as test-particles. We 
adopt a momentum-dependent, Bohm-type diffusion 
coefficient and assume that the magnetic field strength 
scales with the gas density. 

The proton spectrum at the shock, gp{xs,p), and 
the volume-integrated proton spectrum, Gp{p) extends 
to Pmax in equation (9), which increases linearly with 
time. On the other hand, the electron spectrum at 
the shock, ge{xs,p), approaches to the time-asymptotic 
spectrum for the shock age t > t^q in equation (16) . In 
that regime, our time-dependent results with a Bohm- 
type diffusion are qualitatively consistent with the an- 
alytic solutions for a stead-state plane shock, which 
were previously presented by several authors such 
as Heavens & Meisenheimer (1987) with momentum- 
independent K and Zirakashvili & Aharonian (2007) 
with m,om,entum- dependent n{p). So we will re-iterate 
some of the major findings discussed by those authors 
and add new insights obtained from our nonlinear DSA 
simulations. 

1) First of all, we re-derive two characteristic mo- 
menta: the cutoff momentum, poq m equation (15) (for 
the Bohm-type diffusion coefficient) and the break mo- 
mentum, pbr in equation (21). Note that Peq is a time- 
asymptotic quantity that is achieved when the DSA 
energy gain balances the synchrotron/IC energy losses, 
while pbr is a time-dependent quantity that is deter- 
mined by 'aging' of electrons due to synchrotron/IC 
cooling downstream of the shock. 

2) The time-asymptotic electron distribution func- 
tion at the shock, fe(xs,p), has a Gaussian cutoff as 
exp(— p^/Pgq), which agrees well with the analytic form 
suggested by Zirakashvili & Aharonian (2007). 

3) Behind the shock synchrotron/IC cooling domi- 
nates over DSA, so the electron spectrum, fe{x,p), cuts 
off at progressively lower Pcut{d) ~ {u2/DBl2)d~^ , 
which decrease with the distance, d = Xg — x, from the 



shock and is smaller than peq- This cutoff momentum 
is determined by the cooling rate DBI2, independent 
of DSA. 

3) The electron cooling can be represented by the ad- 
vection of the distribution function ge{p) = fe{p)p'^ in 
y = ln(p) space with the advection speed, V = —DB^p 
(see Eq. [3]). This causes the electron spectrum, 
fe{x,p), cuts off more sharply as the distance down- 
stream from the shock increases. 

4) Because the synchrotron/IC cooling time de- 
creases with momentum as frad oc p~^, thickness of 
the electron distribution is inversely proportional to the 
momentum, i.e., Ax{p) = U2 ■ trad oc p^^. Then the 
electron spectrum integrated over to the downstream 
region steepens as -Fe,2(p) cx: for pbr(t) < p < 
Peq, when the spectrum at the shock is fe{xs,p) cx 
The break momentum decreases with the shock age as 
Pbr (xt~^ (see Eq. [21]). 

5) Only highest energy electrons diffuse upstream to 
the distance of d ^ K{pcq) /ug, so the upstream inte- 
grated spcctriim has a much harder spectrum than the 
downstream integrated spectrum and it peaks at Peq- 

6) For a CR modified shock, both proton and elec- 
tron spectra exhibit the well-known concave curva- 
tures. Thus the volume integrated spectrum, F^ip), 
cannot be represented by the canonical broken power- 
law spectrum. In this regime, the radio synchrotron 
index, a, could be steeper than 0.5 even for a high 
sonic Mach number. Also in the case of small Alfvenic 
Mach mmiber {i.e., large Bq and small po), the spec- 
tral slope could be even steeper due to the Alfvenic 
drift effect. Moreover, detail analysis of the X-ray syn- 
chrotron emission near the cutoff frequency may pro- 
vide some information about the effect of nonlinear 
DSA at shocks. Spectral characteristics of the syn- 
chrotron emission from a CR modified shock will be 
presented elsewhere. 
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